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Abstract 

This paper presents a simple, efficient, and high-order accurate sliding-mesh interface approach to the spec¬ 
tral difference (SD) method. We demonstrate the approach by solving the two-dimensional compressible 
Navier-Stokes equations on quadrilateral grids. This approach is an extension of the straight mortar method 
originally designed for stationary domains [7, 8], it employs curved dynamic mortars on sliding-mesh inter¬ 
faces to couple rotating and stationary domains. On the nonconforming sliding-mesh interfaces, the related 
variables are first projected from cell faces to mortars to compute common fluxes, and then the common 
fluxes are projected back from the mortars to the cell faces to ensure conservation. To verify the spatial 
order of accuracy of the sliding-mesh spectral difference (SSD) method, both inviscid and viscous flow cases 
are tested. It is shown that the SSD method preserves the high-order accuracy of the SD method. Mean¬ 
while, the SSD method is found to be very efficient in terms of computational cost. This novel sliding-mesh 
interface method is very suitable for parallel processing with domain decomposition. It can be applied to 
a wide range of problems, such as the aerodynamics of rotorcraft, wind turbines, and oscillating wing wind 
power generators, etc. 

Keywords: spectral difference method, rotating mesh, sliding mesh, high-order methods, unstructured grid 


1. Introduction 

High-order (third and above) numerical methods are becoming more and more popular in recent years due 
to their capability of producing more accurate solutions on relatively coarse grid [31]. The spectral difference 
(SD) method is one discontinuous higlr-order method for solving the conservation laws on unstructured grids 
[15, 32, 28, 11]. This method is an extension of the staggered multi-domain high-order method originally 
designed by Kopriva and Kolias [9]. It was shown that the SD method also has strong connection with 
the Flux Reconstruction/Correction Procedure via Reconstruction (FR/CPR) methods [5], and it shares 
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similarity with the quadrature-free discontinuous Galerkin method [19]. The stability of a particular choice 
of flux points for the SD method was proved by Jameson [6] for the one-dimensional linear wave equation. 
Although the proof has not been generalized to higher-dimensional tensor-product elements, we have not 
observed numerical instability from several successful turbulent flow simulations [14, 1, 22]. 

We have seen more and more applications of the SD method to realistic flow simulations, for example, 
for large eddy simulations on fixed grids [14, 1, 22, 25, 24, 16]. The SD method is also particularly promising 
for simulating vortex-dominated flows on moving and deforming grids [23, 34]. Liang et al. [12] extended 
the SD method for simulating two-dimensional unsteady flows around a plunging or pitching airfoil. DeJong 
and Liang [2] studied three-dimensional vortex induced vibrations using the SD method. 

However, when the mesh undergoes very large rotation motion, such as for flows around rotating pro¬ 
pellers or passing a flapping wing with very large pitching angles, remeshing [29, 30] is required. Our goal is 
to involve the minimum number of remeshing and simultaneously preserve the high-order accuracy of the SD 
method. This motivates us to develop a new approach to the SD method for coupled rotating and stationary 
domains with sliding-mesh interfaces. In our approach, both inviscid and viscous fluxes on the sliding-mesh 
interfaces are constructed using a newly developed curved dynamic mortar method. The mortar method 
on fixed grids was originally proposed for incompressible flows by Mavriplis [18]. Kopriva [7, 8] proved the 
conservation property of the mortar method for the compressible flow equations and applied it to the com¬ 
pressible Euler and Navier-Stokes equations on stationary domains using structured grids. In this paper, 
we show that our sliding-mesli approach is as simple as those designed for low-order numerical methods 
[33, 20] while preserving the high-order accuracy of the SD method. This simple but novel sliding-mesh 
spectral difference (SSD) method can have a wide range of applications, such as rotorcraft aerodynamics, 
wind turbine wake dynamics, and oscillating wing wind power generators. 

The paper is organized as follows: Section 2 gives the two-dimensional compressible Navier-Stokes equa¬ 
tions on stationary and rotating domains. Section 3 reviews the SD method and presents the SSD method 
in detail. Verification studies and applications are reported in Section 4. Section 5 concludes the paper. 


2. The governing equations 

2.1. The compressible Navier-Stokes equations on stationary domain 

We consider the two-dimensional unsteady compressible Navier-Stokes equations in conservative form, 


d>Q dF 

dt dx 
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where Q is the vector of conservative variables, F and G are the x and the y flux vectors. These terms have 
the following expressions, 


Q = [p pu pv E] t , 

( 2 ) 

F = F in V (Q) + F t , is (Q, VQ), 

(3) 

G = G inv (Q) + G v i s (Q,VQ), 

(4) 


where p is the fluid density, u and v are the x and the y velocities, E is the total energy per volume defined 
as E = p/(y — 1) + \p{u 2 + v 2 ), p is the pressure, 7 is the ratio of specific heats and is set to 1.4 (i.e., the 
typical value for the air in standard conditions). 

As shown in Equations (3) and (4), the fluxes have been divided into inviscid and viscous parts. The 
inviscid fluxes are only functions of conservative variables, which are 
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The viscous fluxes are functions of the conservative variables as well as their gradients, 
following expressions, 
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where 77 ,- is the shear stress tensor and is related to the velocity gradients as 77 , = p(uij + Ujp) + A 6ijUk,k, 
p is the dynamic viscosity, A = — 2/3/x based on Stokes’ hypothesis, is the Kronecker delta, k is the 
thermal conductivity, T is the temperature which is related to density and pressure through the ideal gas 
law p = pRT , where R is the gas constant. 


2.2. The compressible Navier-Stokes equations on rotating domain 

On the rotating domains, we implement a simplified equation which is equivalent to the Arbitrary 
Lagrange-Eulerian (ALE) [4] form of Equation (1). Due to grid motion, the inviscid fluxes are modified to 
take the following forms, 
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where u g and v g are the x and the y grid velocities, respectively. The viscous fluxes and all other variables 
are uninfluenced and take the same expressions as those in the previous section. 

For a domain rotating at angular velocity u>, the grid velocities are ( u g ,v g ) = u) x r, where r is the 
position vector with respect to rotating center. For all test cases in the present study, u) is known as a 
priori, thus grid velocities and coordinates are updated analytically on the rotating domains. 


2.3. The transformed equations 

As will be discussed in the next section, we map each quadrilateral grid cell in the physical domain to a 
standard square element in a computational domain. This mapping facilitates the construction of solution 
and flux polynomials. As a result, we only need to solve a set of transformed equations within each standard 
element. Let us assume that the physical coordinates (x,y) are mapped to the computational ones (£, 77 ) 
through a transformation: x = x{f, rj), y = y(£, rf). It can be shown that Equation (1) will take the following 
conservative form after coordinates transformation, 


<9Q <9F <9G 

-I-1- = n 

dt df dr) 

where Q = |J1Q, and the transformed fluxes F, G are related to the physical ones as 
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where \ J\ is determinant of the Jacobian matrix, and J 1 is the inverse Jacobian matrix, 
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3. Numerical methods 

In this section, we first give a brief review of the SD method. Subsequently, we describe a newly 
formulated sliding-mesh interface technique that is built on the SD formulation. For temporal discretization, 
an explicit strong stability preserving Runge-Kutta method [27] is used for all computations throughout this 
paper. 

3.1. The SD method 

For SD method on quadrilateral grids, we first transform each cell in the physical domain to a standard 
square element (0 < £ < 1,0 < 77 < 1) in the computational domain. The transformation can be done 
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through iso-parametric mapping. As was reported by Liang et al. [10, 13], using linear cell defined by four 
nodes is not sufficient for problems involving curved boundaries. Higlr-order cubic cell with twelve nodes 
are used along the curved boundaries to ensure stability and accuracy in the present study. 

After the mapping, solution points (SPs) and flux points (FPs) are defined on each standard element as 
shown in Figure 1 for a third-order scheme. For an N-th order SD method, N SPs are required along each 
coordinate direction to construct degree (N — 1) solution polynomials, and (N + 1) FPs are employed in 
each direction to construct degree N flux polynomials. In the current implementation, the SPs: X s , where 
s = 1,2are chosen as N Chebyshev-Gauss points. The FPs: X s+1 / 2 , where s = 0,1,2are 
chosen as (N — 1) Legendre-Gauss points plus two end points to align in a staggered fashion with the SPs. 

1 


V 


0 

Figure 1: Schematic of the distribution of solution points (circles) and flux points (squares) for a third-order SD scheme. 





ensure conservation. In the current implementation, the Rusanov solver [26] has been used for this purpose. 
The common viscous fluxes are computed from common solution and common gradients, and the detailed 
steps can be found in previous papers by Liang et al. [10, 13]. 

3.2. The sliding-mesh interface approach 

Sliding-mesh interfaces are formed between rotating and stationary meshes. The simplest situation 
involves only one rotating mesh and one stationary mesh as shown in Figure 2. The inner mesh can rotate 
while the outer is fixed, or vice versa. Communication between the stationary and the rotating meshes are 
realized through “mortars”. To make the explanation intuitive, we have scaled the inner mesh in order to 
place mortars in between the two coupled meshes. 



Figure 2: Schematic of the distribution of mortars (hatched) between a rotating and a stationary meshes. 

The mortars are arranged in a counterclockwise order. We refer the inner mesh as left (L) and the outer 
mesh as right (R) with respect to the mortars. To facilitate code implementation and reduce computational 
cost, cell faces on both sides of the sliding-mesh interface have been uniformly meshed. A closer look at 
Figure 2 reveals how mortars and cell faces on the sliding-mesh interface are connected: at each time instant, 
a cell face is connected to two mortars, and each mortar is associated with one left and one right cell faces. 
This cell face and mortar connectivity needs to be updated at every stage of the Runge-Kutta time-stepping 
method. As was discussed in [7] for stationary grid, each cell face can have more than two mortars.Thus, 
our sliding-mesh interface method can also be extended to non-uniform meshes. 

Figure 3 shows a cell face f1 and the attached two mortars Si and S 2 . Each curved mortar is mapped 
to a straight edge 0 < z < 1 through ID iso-parametric mapping. Face fl is mapped to a straight edge 
0 < f < 1 when the associated cell is mapped to a standard square element, thus no extra mapping is 
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required. £ and z are related by 


i = o{t)+s(t)z, 


(17) 


where o(t) is the offset of the mortar relative to the bottom node of fl at time t, and s(t) is the relative 
scaling. For the example shown in Figure 3, we have o\ = 0 and si = / L n for Si, 02 = L =1 / L n and 

S 2 = L“ 2 /L n for S 2 , where L denotes the physical length of the face or the mortar elements. 




“2 




Figure 3: Mapping of curved cell face and mortars to straight ones: left, curved face and mortars in physical domain; right, 
straight face and mortars in computational domain. 


According to Equation (14), the solutions on £1 can be represented as 

N 

Q fl = EQ?^). (!8) 

i=1 

where Qp represents solution at the *-th SP on f l, and hi is the Lagrange basis defined in Equation (12). If 
we define the same set of SPs on 0 < z < 1 for each mortar, then the solutions on each mortar element can 
be reconstructed similarly as 

N 

q e = E^)' ( 1Q ) 

i—1 

where Qf is the solution at the *-th SP on a mortar H. 

The procedure for computing Qf is demonstrated in Figure 4(a). For simplicity, we only show the 
process on the left side of mortar S. To get the solutions, we require that 

[ (Cf’ L (z)-Q n (O)h j (z)dz = 0, i = 1,2, ...,7V. (20) 

Jo 

It was shown in [7] that the above requirement is equivalent to an unweighted L 2 projection. Substituting 
Equations (17)-(19) into the above equation and evaluating it at each SP on H will give a system of linear 
equations. The solution of this system when written in matrix form is 

qH,l = pO^HqQ = m -i s ^HqO i (2!) 


where p 0- *'- is the projection matrix from Q to H, and the elements of the matrices M and S s2 ^“ are 

M id = [ hi(z)hj(z)dz, i,j = 1,2,..., N, (22) 

Jo 

S?j*~ = / hi(o +sz)hj(z)dz, i,j = 1,2,..., N, (23) 

Jo 
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where o and s are the offset and the scaling of H with respect to fi. It is important to note that o and s are 
time-dependent for the sliding-mesh interface method. 

The right solution vector can be computed in the same way. Having both the left and the right 

solutions on a mortar, the Rusanov solver is employed to compute the common inviscid flux F^. This flux 
is then transformed to the computational flux as Ff nv according to Equation (9). 



Figure 4: Projection between face and mortar: (a) from left face to left side of mortar, (b) from two mortars back to the 
associated left face. 


As shown in Figure 4(b), to project the common inviscid fluxes F" n \, and Ff£ v back to face f2, we require 
that, 

r°2 ^ _ cl ^ ^ 

/ (F« w (0-F % v (z))h J (Qdt+ F? nv (O--F% v (z))h j (QdZ = 0, j = 1,2,..., N, (24) 

JO J 02 

where F^(£) is the inviscid flux polynomial on face H. The solution of the above equation when written 
in matrix form is 


_ p- 

inv 


+ P" 


2 F“?„ = 8 1 M~ 1 S S 


s 2 M~ 1 S' 


(25) 


where the matrix M is identical to that of Equation (21), the matrices and S“ 2_>f2 are simply the 

transposes of and respectively. 

For the computation of the common viscous fluxes, we first compute the common solution on each mortar 
as the average of the left and the right solutions, 

Q S = l(Q S ’ L + Q S ’ R )- ( 26 ) 

This common solution is then projected back to the cell faces in the same procedure as for the inviscid flux 
in Equation (25). After that, solution gradients and viscous fluxes are updated on the cell faces on both 
sides of the interface. The viscous fluxes FR s on the cell faces are projected to the mortars in the same way 
as Equation (21). The common viscous flux Fj; is on a mortar is taken as the average of the left and the 
right viscous fluxes, 

pH _ l/pH.L pS,R\ ,< 27 ) 

VIS 2 ' VIS 1 VIS /• V 1 / 

The final step is to project F“ is back to cell faces, which is identical to the process showed in Equation (25). 



Since a uniform mesh is used for the cell faces on the sliding-mesh interface, the matrix S only needs 
to be computed for the first two mortars, and can be reused by other corresponding mortars. At the same 
time since the matrix M is time independent, it can be computed and stored before the actual calculation. 
To compute the integrals in Equations (22) and (23), one can use the Clenshaw-Curtis quadrature method 
as was used in [7]. In this study, the integrand is casted into a general form as a product of 2 (N — 1) first 
degree polynomials, and we implement a recursive algorithm to compute the integrals analytically. This 
approach requires the least number of operations which is much more efficient than the Clenshaw-Curtis 
quadrature method. 


4. Numerical tests 

In this section we test the spatial accuracy of the SSD method on both inviscicl and viscous flows, and 
then apply this method to study an external and an internal flows. A five-stage fourth-order Runge-Kutta 
method for time stepping [27] is used for all test cases. In each test case for demonstrating the orders of 
accuracy of spatial discretizations, the time step size is reduced successively until the final errors do not 
change with it. This ensures that the temporal discretization errors are negligible and the final errors can 
represent the spatial discretization errors. 


4-1. Euler vortex flow 

Euler vortex flow has been widely used to test the accuracies of inviscicl flow solvers [3, 31]. In this 
problem, an isentropic vortex is superimposed to and convected by a uniform mean flow. The Euler vortex 
flow in an infinite domain at time t can be analytically described as 
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where U 00 , p 00 , Poo , M a0 are the mean flow speed, density, pressure and Mach number, respectively. 9 is 
the direction of the mean flow (i.e. the direction along which the vortex is convected), e and r c can be 
interpreted as the vortex strength and size. The relative coordinates (. x r ,y r ) are defined as 


x r = x — Xo — ut, 
yr = y-yo- vt, 
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(32) 

(33) 



where u = Uao cos 9, v = Uao sin 8 are the x and y components of the mean velocity, (xo,yo) is the initial 
position of the vortex. The analytical solution of the Euler vortex flow problem within a square domain 
(0 < x, y < L) with periodic boundary conditions can be achieved by replacing the relative coordinates with 
the following expressions, 

x r =x r - [—-— J ■ L, (34) 

Vr=Vr- L-£-J ' L) ( 35 ) 

where the floor operator |_£j gives the largest integer that is not greater than a real number x. The x r and 
y r on the right hand sides of Equations (34) and (35) are from Equations (32) and (33). 

In this test, the uniform mean flow is chosen as (Uao, Poo,Poo) = (1, 1, 1) with a Mach number of = 
0.3. The flow direction is set to 9 = arctan(l/2). A vortex with parameters: e = 1, r c = 1, is superimposed 
to the mean flow. The domain size is 0 < x,y < 10 (i.e. L = 10), and the vortex is initially located at the 
domain center ( xo,yo ) = (5,5). Periodic boundary conditions are applied in both x and y directions. 

Figure (5) shows a computational mesh with 700 cells. The mesh has been decomposed into two parts: a 
rotating inner part with a radius of 2; a fixed outer part which takes the rest of the computational domain. 
Three meshes with 180, 700 and 2731 cells have been used for accuracy tests. For all three cases, the inner 
part is set to rotate at an angular speed of to = 1.0. 



Figure 5: Mesh with 700 cells at a time instant for the Euler vortex flow simulation (blue circle indicates sliding-mesh interface). 

Figure (6) compares the density contours obtained with the fourth-oder SSD method on the finest mesh 
with the exact solution at t = 2. As we can see, the solver resolves the vortex very well, and we see almost 
no visible difference between the exact solution and the numerical one. 
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Figure 6: Contours of density at the time instant t — 2 for the Euler vortex flow. Left, exact solution; right, numerical solution 
from fourth-order scheme (blue circle indicates location of sliding interface). 

Furthermore, Tables 1 and 2 give the spatial accuracy of the scheme, where the LI and L2 errors are 
computed from density at t = 2 when vortex center is traveled right onto the sliding-mesh interface. From 
the two tables we see that the SSD method gives very reasonable order of accuracy. 


cells 

LI error 

order 

L2 error 

order 

180 

3.16E-4 

- 

8.13E-4 

- 

700 

4.58E-5 

2.85 

1.13E-4 

2.90 

2731 

7.00E-6 

2.80 

1.73E-5 

2.83 


Table 1: Errors and orders of accuracy of the third-order scheme for the Euler vortex flow simulation. 


cells 

LI error 

order 

L2 error 

order 

180 

5.43E-5 

- 

1.26E-4 

- 

700 

3.27E-6 

4.14 

8.02E-6 

4.06 

2731 

2.26E-7 

4.03 

5.50E-7 

4.00 


Table 2: Errors and orders of accuracy of the fourth-order scheme for the Euler vortex flow simulation. 


To see how efficient the SSD method is, we compare the total computational time and the communication 
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time on the sliding-mesh interface in Tables 3 and 4 for the third- and fourth-order schemes, respectively. 
Times in both tables are collected for 100 computational steps and do not include any post-processing time. 
It is seen that for all test cases, communication on the sliding-mesh interface takes only a few percent of the 
total computational time, which clearly shows that the SSD method is efficient. It is interesting that the 
relative communication time (represented by the percentage) decreases as either the number of cells or the 
order of schemes increases. This is due to the fact that cells are one dimension higher than faces: when we 
perform a mesh refinement, the total number of faces in the domain grows faster than on the sliding-mesh 
interface; when we increase the scheme order, the total number of degrees of freedom in the domain also 
grows faster than on the sliding-mesh interface. 


cells 

total time 

comm, time 

percentage 

180 

0.254944 

0.017657 

6.92% 

700 

1.013039 

0.041758 

4.12% 

2731 

4.607205 

0.097717 

2 .12% 


Tabic 3: Total computation time and interface communication time (both in seconds) for 100 computational steps using a 
third-order scheme for the Euler vortex flow simulation. 


cells 

total time 

comm, time 

percentage 

180 

0.420763 

0.023282 

5.53% 

700 

1.763656 

0.058833 

3.34% 

2731 

7.511551 

0.125820 

1 .68% 


Table 4: Total computation time and interface communication time (both in seconds) for 100 computational steps using a 
fourth-order scheme for the Euler vortex flow simulation. 


4-2. Taylor- Couette flow 

To test the order accuracy on viscous flow, we use the laminar Taylor-Couette flow as the test case. 
Previous researchers such as Liang et al. [13], Michalak and Ollivier-Gooch [21] used similar flows to test 
the accuracy of their solvers. In the present test, the inner cylinder has a radius of r, = 1, the outer cylinder 
has radius of r 0 = 2. Both boundaries are set to be isothermal walls. The domain has been divided into two 
parts at r = 1.5. The inner part rotates at an angular speed of w, = 1 while the outer part stays stationary. 
The Reynolds number based the inner cylinder radius and speed is Re = 10. The Mach number on the inner 
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wall is set to M t = 0.1. Three meshes with 192, 768 and 3072 cells are used for the tests. Figure 7 shows 
the mesh with 192 cells. 

Figure 8 shows the steady state contours of the u velocity component and the Mach number obtained 
with a fourth-oder scheme on the finest mesh. We see that the Mach contours at the steady state are a 
series of concentric circles, and the u velocity is highly symmetric. These results are consistent with our 
expectations. 

The exact solution for the circumferential velocity has the following relation to radius r, 

r 0 /r - r/r 0 , . 

ve=UiTi — - -—, (36) 

r 0 /ri - ri/r 0 

The x component of this velocity (i.e., u ) is used to compute the LI and L2 error norms. From Table 5 and 
Table 6 we see that the SSD method preserves the high-order accuracy for viscous flow as well. 



Figure 7: Mesh with 192 cells at a time instant for the Taylor-Couette flow simulation (blue circle indicates sliding-mesh 
interface). 


cells 

LI error 

order 

L2 error 

order 

192 

5.90E-5 

- 

8.71E-5 

- 

768 

6.95E-6 

3.09 

9.82E-6 

3.15 

3072 

8.43E-7 

3.07 

1.09E-6 

3.16 


Table 5: Errors and orders of accuracy of a third-order scheme for the Taylor-Couette flow simulation. 


The total computational time and the sliding-mesh interface communication time are shown in Table 7 
and Table 8 for third- and fourth-order schemes, respectively. Again, data in both tables is collected for 100 
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Figure 8: Contours of u velocity component and the Mach number for the Taylor-Couette flow (dashed circle indicates location 
of sliding-mesh interface). 


cells 

LI error 

order 

L2 error 

order 

192 

9.72E-6 

- 

1.52E-5 

- 

768 

6.16E-7 

3.98 

1.01E-6 

3.91 

3072 

4.22E-8 

3.92 

6.60E-8 

3.92 


Table 6: Errors and orders of accuracy of a fourth-order scheme for the Taylor-Couette flow simulation. 

computational steps and do not include any post-processing time. It is seen that for viscous flow simulations 
the SSD method remains very efficient. 


cells 

total time 

comm, time 

percentage 

192 

1.267700 

0.105698 

8.34% 

768 

4.394639 

0.234498 

5.34% 

3072 

18.138023 

0.561518 

3.10% 


Tabic 7: Total computation time and interface communication time (both in seconds) for 100 computational steps using a 
third-order scheme for the Taylor-Couette flow simulation. 
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cells 

total time 

comm, time 

percentage 

192 

2.075053 

0.139513 

6.92% 

768 

7.628162 

0.322616 

4.12% 

3072 

30.852314 

0.729768 

2.12% 


Table 8: Total computation time and interface communication time (both in seconds) for 100 computational steps using a 
fourth-order scheme for the Taylor-Couette flow simulation. 


4-3. Flow over a rotating elliptic cylinder 

To further verify the approach for flow involving complex geometries, we simulate a laminar flow over a 
two-dimensional elliptic cylinder in this section. Maruoka [17] and Zhang et al. [35] studied incompressible 
flow over a rotating elliptic cylinder using the finite element and finite volume methods, respectively. Both 
studies use Chimera grids for communication between the foreground rotating mesh and the background 
stationary mesh. The freestream Mach number in the present test is set to = 0.05 to mitigate the 
compressibility effects in order to compare with their incompressible flow results. 

The elliptic cylinder has a major and a minor axis lengths of 1.0 and 0.5, respectively. Initially, the major 
axis is parallel to the freestream. The cylinder rotates counterclockwisely at an angular speed of c o = 0.5-7T. 
The Reynolds number based on the freestream velocity and the major axis length is Re = 200. Figure 
9 shows a schematic of the computational domain. Slip boundary conditions are applied on the top and 
bottom boundaries. Dirichlet boundary condition is used at the inlet, and fixed pressure boundary condition 
is applied at the outlet. Finally, no-slip isothermal wall boundary condition is used on the cylinder surface. 



Figure 9: Schematic of the computational domain for flow over a rotating elliptic cylinder (not to scale). 
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The inner rotating domain has a radius of 1.5 and is meshed with 1280 cells. The rest of the domain is 
stationary and has 7391 cells. Mesh refinement are performed around the leading and the trailing edges, as 
well as in the wake region. Figure 10 shows two local views of the mesh. The nondimensional time step size 
AtUoo/L for the simulation is set to 1.0 x 10~ 4 , where L is the major axis length and Uoo is the freestream 
velocity. 



Figure 10: Two local views of the mesh for simulation of flow around a rotating elliptic cylinder (blue circle indicates sliding- 
mesh interface). 

Both third- and fourth-order schemes were tested for this flow and no visible difference was observed 
between the solutions. For this reason, we only present results from the fourth-order scheme. As was 
noticed by Maruoka [17] and Zhang et al. [35], the fully developed flow takes a periodic pattern as the 
cylinder rotates. The lift and the drag coefficients in one rotating period are shown in Figure 11. It is seen 
that the present results agree very well with the previously published results. 

Figure 12 shows the streamlines superimposed on vorticity contours at a series of time instants in one 
rotating period. A clockwise and a counterclockwise vortices appear alternatively around the two ends of 
the cylinder. From Figure 12 (h) and Figure 12 (a), we see that a clockwise vortex is formed at the cylinder 
leading edge as the cylinder rotates, and this vortex is then shed off from the leading edge and travels 
downstream to hit the trailing edge. From Figure 12 (b) to Figure 12 (d), the same clockwise vortex is again 
shed off from the trailing edge and then convected towards downstream. From Figure 12 (e) to Figure 12 
(g), a counterclockwise vortex slowly emerges around the other end of the cylinder and is then convected 
downstream without reattaching to the cylinder. This process repeats as the cylinder rotates, and a vortex 
street forms downstream of the cylinder. 

The efficiency of the SSD method is shown in Table 9 for this case. The results in the table confirm 
our previous conclusion, i.e., the relative communication time generally decreases as the number of cells or 
the order of scheme increases. In fact, the interface communication time in the table is almost negligible 
comparing to the total computational time. 
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Figure 11: Lift and drag coefficients for flow over an elliptic cylinder. 


order 

total time 

comm, time 

percentage 

3 

38.924785 

0.143083 

0.37% 

4 

70.655175 

0.185890 

0.26% 


Table 9: Total computation time and interface communication time (both in seconds) for 100 computational steps for simulation 
of flow over a rotating elliptic cylinder. 

4-4■ Flow inside a 2D stirred tank 

In this last test case, we apply the SSD method to simulate laminar flow inside a 2D stirred tank. 
Figure 13 shows the unstructured quadrilateral mesh used for this case. As we can see, the tank has several 
components: an inner circular wall with a radius of 0.5; an outer circular wall with a radius of 5; six 
uniformly distributed agitating blades, each extends from r = 1 to r = 2 and has a thickness of 0.1; four 
baffles installed on the outer wall, each of them has a height of 1 and a thickness of 0.1. The computational 
domain is split into an inner rotating part and an outer fixed part, resulting in a sliding-mesh interface at 
r = 3. Mesh has been refined around the blades, baffles, and the wall boundaries. The resulted mesh has 
14,990 cells in total. 

In this simulation, the inner circular wall and the blades rotate at an angular speed of to = 1. The 
Reynolds number base on the inner wall diameter and the angular speed is Re = 100. Flow on the inner 
wall surface has a Mach number of Mj = 0.1. No slip isothermal wall boundary conditions are used on the 
inner wall surface, the baffles, and the outer walls. Adiabatic wall boundary conditions are used on the six 
blades. A time step size of At = 1.0 x 10~ 4 is used. 
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Figure 12: Streamlines and vorticity contours (blue color means negative value, red color means positive) for flow over a rotating 
elliptic cylinder. 


Initially, the flow field is set to be uniform and stationary. Figure 14 shows the flow fields at four different 
states by visualizing the density. It is seen in Figure 14(a) that at the initial transient state, the flow behaves 
like flow around a propeller: fluid is “squeezed” and “pushed” away from the blades, resulting in a concentric 
flow pattern, and large fluctuations are observed in the flow field. Quickly after the flow reaches the baffles 
and the outer wall, the flow field becomes very chaotic: vortical structures generated by flow passing the 
baffles, bouncing pressure waves, blades induced vortices, unsteady boundary layers, etc., as can be seen 
in Figure 14(b). As the blades continue rotating for a longer time, the chaotic flow structures are slowly 
being dissipated, and organized large flow structures emerge as shown in Figure 14(c). Finally the flow field 
reaches a quasi steady state in the rotating reference frame, carrying little variation with time, as shown 
in Figure 14(d). Interestingly, the flow is becoming more and more uniform at the later stage of stirring. 
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Figure 13: Mesh for simulation of flow inside a 2D stirred tank (blue circle indicates sliding-mesh interface). 

Variations on the density field get smaller and smaller when comparing Figure 14 (a)-(d) along time. 

As for all previous cases, we monitored the efficiency of the SSD method for this case. The results are 
shown in Table 10. Again, as can be seen from the table, the SSD method is very efficient with negligible 
communication time. 


order 

total time 

comm, time 

percentage 

3 

85.547403 

0.328987 

0.38% 

4 

148.780491 

0.437902 

0.29% 


Table 10: Total computation time and interface communication time (both in seconds) for 100 computational steps for simu¬ 
lation of flow inside a 2D stirred tank. 


5. Conclusions 

In this paper, a novel, simple, efficient, and high-order accurate sliding-mesh interface method is reported 
for subsonic compressible flows. The sliding-mesh spectral difference (SSD) method has been successfully 
developed and tested for several inviscid and viscous flow problems. The SSD method retains the high-order 
accuracy of the SD method. It is also shown that the SSD method is very efficient as it introduces negligible 
extra computational cost to the SD method for realistic flow simulations. In this paper, we demonstrated 
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Figure 14: Density contours of flow inside a 2D stirred tank at four distinct states. 


the approach on uniformly meshed sliding-mesh interfaces where each cell face has two mortars. The sliding- 
mesh interface method can be extended to handle more than two mortars for each face. The SSD method is 
also very suitable for parallel computing. Since there is no overlapping in grids and the sliding-mesh interface 
method introduces negligible computational time, thus each domain can be decomposed and distributed to 
the processors to achieve load balancing. Finally, this high-order curved sliding-mesh interface method can 
also be extended to other discontinuous high-order methods for compressible flows. 
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